home *** CD-ROM | disk | FTP | other *** search
/ Sprite 1984 - 1993 / Sprite 1984 - 1993.iso / src / lib / m / j0.c < prev    next >
C/C++ Source or Header  |  1988-07-11  |  5KB  |  204 lines
  1. /*
  2.  * Copyright (c) 1985 Regents of the University of California.
  3.  * All rights reserved.  The Berkeley software License Agreement
  4.  * specifies the terms and conditions for redistribution.
  5.  */
  6.  
  7. #ifndef lint
  8. static char sccsid[] = "@(#)j0.c    5.2 (Berkeley) 4/29/88";
  9. #endif /* not lint */
  10.  
  11. /*
  12.     floating point Bessel's function
  13.     of the first and second kinds
  14.     of order zero
  15.  
  16.     j0(x) returns the value of J0(x)
  17.     for all real values of x.
  18.  
  19.     There are no error returns.
  20.     Calls sin, cos, sqrt.
  21.  
  22.     There is a niggling bug in J0 which
  23.     causes errors up to 2e-16 for x in the
  24.     interval [-8,8].
  25.     The bug is caused by an inappropriate order
  26.     of summation of the series.  rhm will fix it
  27.     someday.
  28.  
  29.     Coefficients are from Hart & Cheney.
  30.     #5849 (19.22D)
  31.     #6549 (19.25D)
  32.     #6949 (19.41D)
  33.  
  34.     y0(x) returns the value of Y0(x)
  35.     for positive real values of x.
  36.     For x<=0, if on the VAX, error number EDOM is set and
  37.     the reserved operand fault is generated;
  38.     otherwise (an IEEE machine) an invalid operation is performed.
  39.  
  40.     Calls sin, cos, sqrt, log, j0.
  41.  
  42.     The values of Y0 have not been checked
  43.     to more than ten places.
  44.  
  45.     Coefficients are from Hart & Cheney.
  46.     #6245 (18.78D)
  47.     #6549 (19.25D)
  48.     #6949 (19.41D)
  49. */
  50.  
  51. #include <math.h>
  52. #if defined(vax)||defined(tahoe)
  53. #include <errno.h>
  54. #else    /* defined(vax)||defined(tahoe) */
  55. static double zero = 0.e0;
  56. #endif    /* defined(vax)||defined(tahoe) */
  57. static double pzero, qzero;
  58. static double tpi    = .6366197723675813430755350535e0;
  59. static double pio4    = .7853981633974483096156608458e0;
  60. static double p1[] = {
  61.     0.4933787251794133561816813446e21,
  62.     -.1179157629107610536038440800e21,
  63.     0.6382059341072356562289432465e19,
  64.     -.1367620353088171386865416609e18,
  65.     0.1434354939140344111664316553e16,
  66.     -.8085222034853793871199468171e13,
  67.     0.2507158285536881945555156435e11,
  68.     -.4050412371833132706360663322e8,
  69.     0.2685786856980014981415848441e5,
  70. };
  71. static double q1[] = {
  72.     0.4933787251794133562113278438e21,
  73.     0.5428918384092285160200195092e19,
  74.     0.3024635616709462698627330784e17,
  75.     0.1127756739679798507056031594e15,
  76.     0.3123043114941213172572469442e12,
  77.     0.6699987672982239671814028660e9,
  78.     0.1114636098462985378182402543e7,
  79.     0.1363063652328970604442810507e4,
  80.     1.0
  81. };
  82. static double p2[] = {
  83.     0.5393485083869438325262122897e7,
  84.     0.1233238476817638145232406055e8,
  85.     0.8413041456550439208464315611e7,
  86.     0.2016135283049983642487182349e7,
  87.     0.1539826532623911470917825993e6,
  88.     0.2485271928957404011288128951e4,
  89.     0.0,
  90. };
  91. static double q2[] = {
  92.     0.5393485083869438325560444960e7,
  93.     0.1233831022786324960844856182e8,
  94.     0.8426449050629797331554404810e7,
  95.     0.2025066801570134013891035236e7,
  96.     0.1560017276940030940592769933e6,
  97.     0.2615700736920839685159081813e4,
  98.     1.0,
  99. };
  100. static double p3[] = {
  101.     -.3984617357595222463506790588e4,
  102.     -.1038141698748464093880530341e5,
  103.     -.8239066313485606568803548860e4,
  104.     -.2365956170779108192723612816e4,
  105.     -.2262630641933704113967255053e3,
  106.     -.4887199395841261531199129300e1,
  107.     0.0,
  108. };
  109. static double q3[] = {
  110.     0.2550155108860942382983170882e6,
  111.     0.6667454239319826986004038103e6,
  112.     0.5332913634216897168722255057e6,
  113.     0.1560213206679291652539287109e6,
  114.     0.1570489191515395519392882766e5,
  115.     0.4087714673983499223402830260e3,
  116.     1.0,
  117. };
  118. static double p4[] = {
  119.     -.2750286678629109583701933175e20,
  120.     0.6587473275719554925999402049e20,
  121.     -.5247065581112764941297350814e19,
  122.     0.1375624316399344078571335453e18,
  123.     -.1648605817185729473122082537e16,
  124.     0.1025520859686394284509167421e14,
  125.     -.3436371222979040378171030138e11,
  126.     0.5915213465686889654273830069e8,
  127.     -.4137035497933148554125235152e5,
  128. };
  129. static double q4[] = {
  130.     0.3726458838986165881989980e21,
  131.     0.4192417043410839973904769661e19,
  132.     0.2392883043499781857439356652e17,
  133.     0.9162038034075185262489147968e14,
  134.     0.2613065755041081249568482092e12,
  135.     0.5795122640700729537480087915e9,
  136.     0.1001702641288906265666651753e7,
  137.     0.1282452772478993804176329391e4,
  138.     1.0,
  139. };
  140.  
  141. double
  142. j0(arg) double arg;{
  143.     double argsq, n, d;
  144.     double sin(), cos(), sqrt();
  145.     int i;
  146.  
  147.     if(arg < 0.) arg = -arg;
  148.     if(arg > 8.){
  149.         asympt(arg);
  150.         n = arg - pio4;
  151.         return(sqrt(tpi/arg)*(pzero*cos(n) - qzero*sin(n)));
  152.     }
  153.     argsq = arg*arg;
  154.     for(n=0,d=0,i=8;i>=0;i--){
  155.         n = n*argsq + p1[i];
  156.         d = d*argsq + q1[i];
  157.     }
  158.     return(n/d);
  159. }
  160.  
  161. double
  162. y0(arg) double arg;{
  163.     double argsq, n, d;
  164.     double sin(), cos(), sqrt(), log(), j0();
  165.     int i;
  166.  
  167.     if(arg <= 0.){
  168. #if defined(vax)||defined(tahoe)
  169.         extern double infnan();
  170.         return(infnan(EDOM));        /* NaN */
  171. #else    /* defined(vax)||defined(tahoe) */
  172.         return(zero/zero);    /* IEEE machines: invalid operation */
  173. #endif    /* defined(vax)||defined(tahoe) */
  174.     }
  175.     if(arg > 8.){
  176.         asympt(arg);
  177.         n = arg - pio4;
  178.         return(sqrt(tpi/arg)*(pzero*sin(n) + qzero*cos(n)));
  179.     }
  180.     argsq = arg*arg;
  181.     for(n=0,d=0,i=8;i>=0;i--){
  182.         n = n*argsq + p4[i];
  183.         d = d*argsq + q4[i];
  184.     }
  185.     return(n/d + tpi*j0(arg)*log(arg));
  186. }
  187.  
  188. static
  189. asympt(arg) double arg;{
  190.     double zsq, n, d;
  191.     int i;
  192.     zsq = 64./(arg*arg);
  193.     for(n=0,d=0,i=6;i>=0;i--){
  194.         n = n*zsq + p2[i];
  195.         d = d*zsq + q2[i];
  196.     }
  197.     pzero = n/d;
  198.     for(n=0,d=0,i=6;i>=0;i--){
  199.         n = n*zsq + p3[i];
  200.         d = d*zsq + q3[i];
  201.     }
  202.     qzero = (8./arg)*(n/d);
  203. }
  204.